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Relative frequencies of mechanically stable (MS) packings of frictionless bidisperse disks are stud- 
ied numerically in small systems. The packings are created by successively compressing or de- 
compressing a system of soft purely repulsive disks, followed by energy minimization, until only 
infinitesimal particle overlaps remain. For systems of up to 14 particles most of the MS packings 
were generated. We find that the packings are not equally probable as has been assumed in recent 
thermodynamic descriptions of granular systems. Instead, the frequency distribution, averaged over 
each packing-fraction interval A</), grows exponentially with increasing 0. Moreover, within each 
packing-fraction interval MS packings occur with frequencies that differ by many orders of magni- 
tude. Also, key features of the frequency distribution do not change when we significantly alter the 
packing-generation algorithm — for example frequent packings remain frequent and rare ones remain 
rare. 

These results indicate that the frequency distribution of MS packings is strongly infiuenced by 
geometrical properties of the multidimensional configuration space. By adding thermal fiuctuations 
to a set of the MS packings, we were able to examine a number of local features of configuration 
space near each packing. We measured the time required for a given packing to break to a distinct 
one, which enabled us to estimate the energy barriers that separate one packing from another. We 
found a gross positive correlation between the packing frequencies and the heights of the lowest 
energy barriers eo; however, there is significant scatter in the data. We also examined displacement 
fluctuations away from the MS packings to assess the size and shape of the local basins near each 
packing. The displacement modes scale as di ^ Eq' with 7; ranging from approximately 0.6 for the 
largest eigenvalues to 1.0 for the smallest ones. These scalings suggest that the packing frequencies 
are not determined by the local volume of configuration space near each packing, which would 
require that the dependence of on eo is much stronger than the dependence we observe. The 
scatter in our data implies that in addition to eo there are also other, as yet undetermined variables 
that influence the packing probabilities. 

PACS numbers: 81.05.Rm, 81.05.Kf 83.80.Fg 
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I. INTRODUCTION 

Despite intense study over the past several decades, 
glassy and amorphous materials are stiU poorly under- 
stood. For example, a fundamental explanation for the 
stupendous rise in the viscosity of fragile glass-forming 
liquids as the temperature is lowered near the glass tran- 
sition is still lacking 0. Also, the response of glassy 
and amorphous systems to applied stress is difficult to 
predict because these systems display complex spatio- 
temporal dynamics, such as shear bands strongly 
non-affine and cooperative motion |^ , and dynamical 
heterogeneities 0|. Even basic questions concerning the 
nature of stress and structural relaxation have not been 
adequately addressed. Important open questions include 
1) what are the characteristic rearrangement events that 
lead to stress and structural relaxation, 2) how many 
particles are involved in such rearrangement events, and 
3) are these events correlated and over what length and 
time scales? 

An extremely useful concept for understanding the dy- 
namical and mechanical properties of glassy systems has 
been the potential energy landscape (PEL) formalism f^. 



The PEL is the highly multi-dimensional potential en- 
ergy function that depends on all of the configurational 
degrees of freedom of the system. It has been shown 
that the eq uation of state 0, and dynamical quanti- 
ties 0, 01 such as the diffusion constant and viscosity 
of supercooled liquids and glasses can be calculated in 
terms of geometrical features of the PEL like local poten- 
tial energy minima and low-order saddle points. Related 
studies have also been performed on hard spheres to un- 
derstand the glass transition in these model systems. A 
significant focus of this research has been to explain ki- 
netic arrest in terms of decreasing free volume |l2l 
and configurational entropy [T^ ITsjl near the glass tran- 
sition. 

In this article we investigate particle packings and fea- 
tures of the potential energy landscape in 2d bidisperse 
systems of frictionless disks. The disks interact via a 
finite-range continuous repulsive potential. We focus on 
mechanically stable (MS) packings with vanishing parti- 
cle overlaps. In these mechanically stable packings any 
particle displacement results in an increase of the poten- 
tial energy, i.e., leads to overlap between particles. Thus 
the set of MS packings in our system is equivalent to 
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the set of collectively jammed states |16| for hard disks. 
However, since we consider particles that interact via a 
continuous potential, we can explore not only geometri- 
cal features of configuration space in the neighborhood 
of a given collectively jammed state, but also properties 
of the energy landscape near a given mechanically stable 
packing. 

An important feature of the present work (and our 
related earlier study is that we focus on small sys- 
tems containing 14 or fewer particles. Related studies of 
small hard disk systems have been carried out previously 
[Tsl [Tflj l. but these did not address questions concern- 
ing the frequency with which MS packings occur, which 
is the main focus of this work. We confine our studies 
to small systems for several reasons. First, we believe 
that understanding properties of small nearly jammed 
systems is crucial to developing a theoretical explanation 
for slow stress and structural relaxation in large glassy 
and amorphous systems. Second, in small systems we 
are able to generate nearly all of the mechanically stable 
disk packings, which enables us to accurately measure the 
frequency with which different MS packings occur. Also, 
detailed results obtained for small systems of different 
size can be extrapolated and used to infer properties of 
glassy materials in the large-system limit. 

Results from a number of recent studies emphasize 
that small subsystems are important for understanding 
the slow dynamics displayed in supercooled liquids and 
glasses. For example, experiments on colloidal glasses 
|20l l2l| show that slow relaxation in glassy materials 
occurs through local cage breaking events. Caging be- 
havior has also been shown in a number of computer 
simulations of hard and soft particles, for example in 
Refs. miillll. 

Another indication that small subsystems play an im- 
portant role in determining the dynamics of glass- forming 
liquids can be found in the results of recent numerical 
simulations of large binary hard-disk systems in Ref. (2^ . 
Using an appropriately defined signature of the kinetic 
glass transition, these authors have shown that a slowly 
quenched system falls out of equilibrium at packing frac- 
tions (j) in the range 0.77 < < 0.8 (the larger packing 
fractions correspond to slower quench rates). As revealed 
by our recent study (see also the results presented 
herein), this is the packing- fraction range where the max- 
imum of the density of MS packings occurs for small sys- 
tems of > 12 particles. We believe that this is not a 
coincidence. 

The important role of small subsystems in dense hard- 
particle fluids stems from the fact that small systems 
become geometrically jammed (i.e., the system cannot 
rearrange) at packing fractions where significant free vol- 
ume is still available for the motion of particles in their 
local cages. In contrast, for macroscopic systems even 
infinitesimal free volume per particle makes structural 
rearrangements geometrically possible. 

Such observations suggest that one should seek macro- 
scopic descriptions of glass-forming materials in terms of 



many coupled, nearly jammed small subdomains. At a 
packing fraction at which small subdomains become ge- 
ometrically jammed, rearrangements can occur only if a 
domain is sufficiently uncompressed. Large local density 
fluctuations, however, occur infrequently, because of the 
low compressibility of the surrounding domains at pack- 
ing fractions above the dynamic glass transition. Hence, 
the structural relaxation time is extremely large. A pic- 
ture similar to the one put forward by Adam and Gibbs 
over 40 years ago can also be constructed for systems 
of particles interacting via continuous potentials. In this 
case, relevant packing fractions can be defined using a 
temperature-dependent effective particle diameter. 

A small subsystem of a macroscopic system corre- 
sponds to a subspace in configuration space. Thus, an 
approach that describes the dynamics of glass-forming 
liquids in terms of a set of small coupled subsystems 
should help to reconcile the local picture of glassy mate- 
rials (e.g., trap models j27,]) with the global PEL view. 
However, before predictive theories can be constructed 
based on the ideas outlined above, one needs to obtain 
a more complete understanding of the behavior of small 
nearly jammed systems. 

An analysis of MS packings in small systems is rel- 
evant not only for understanding the slow dynamics of 
glass-forming liquids, but also for explaining the meaning 
of random-close packing in macroscopic athermal amor- 
phous materials. We now turn to this aspect of the prob- 
lem. 

In MS packings composed of touching (but not overlap- 
ping) disks, the number of degrees of freedom is less than 
or equal to the number of constraints. Thus, MS pack- 
ings, or collectively jammed states, can be represented 
by single points in configuration space. In small systems 
we are able to generate nearly all mechanically stable 
disk packings. Thus, we can investigate, separately, the 
number of distinct packings that exist in a given interval 
of packing fraction (i.e., the density of MS packings) 
and the probability with which these packings occur for 
a given generation protocol 0| . 

An analysis of the protocol-independent density of MS 
packings and protocol-dependent probabilities for each 
distinct packing allows us to address the question of why 
random close packing seems to be a well defined quan- 
tity (i.e., many distinct generation protocols give similar 
values for it, for example in Refs. j2^, '2^, "s^, ^T] ) , but it 
is also a poorly defined concept because other protocols, 
for example, those that allow slow thermal quench rates, 
yield a continuous range of packing fractions at which 
jammed states occur 32]. 

In our opinion, questions concerning random close 
packing have not yet been fully resolved. In particular, 
although the maximally random jammed (MRJ) state as 
described in Ref. [s^l is a useful and important concept, 
we believe that it is not the final answer. For example, 
MRJ involves an arbitrary choice of how to characterize 
order in the system and does not address the question 
of why a wide class of generation protocols yield similar 
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(a) (b) 

FIG. 1: Snapshot of (a) a rare mechanically stable packing 
at « 0.82 and (b) another MS packing cf) ^ 0.83 that is 10*^ 
times more frequent. 

values for random-close packing. 

In our recent paper mT] we argued that the random 
close packing density corresponds to the position of the 
peak in the density of MS packings p{(j)) , which becomes a 
(5- function in the infinite system-size limit. As opposed to 
the protocol-dependent probability density P{(l>) for ob- 
taining a packing at a given cj), the density of MS pack- 
ings is a protocol-independent quantity. Unless a given 
protocol is strongly biased toward states in the tail of 
the distribution p{(j)), for example those protocols that 
involve significant thermalization, a well-defined random 
close packed density is obtained in the large system limit. 
This explains why many quite distinct protocols consis- 
tently yield similar values for the random close packed 
density 0rcp- 

This explanation, while quite plausible, leaves a num- 
ber of important issues unresolved. For example, we find 
that individual MS packings occur with extremely dif- 
ferent frequencies even for a typical fast-quench protocol 
|l7j | . Two MS packings at approximately the same pack- 
ing fraction (j) but with drastically different frequencies 
are shown in Fig.^ There are no striking structural dif- 
ferences between these two MS packings, yet, the packing 
shown in Fig.^a) is 10^ times less frequent than that in 
Fig. n^b). Moreover, we find that there are many more 
infrequent MS packings than frequent ones. 

An important issue is what determines whether a par- 
ticular MS packing is frequent or rare. Although it is 
clear that the particular protocol chosen to generate the 
MS packings plays a role in determining the frequency 
distribution p33 |. we argue that geometrical features of 
the PEL also strongly influence the frequency distribu- 
tion. First, our protocols, which involve a sequence of 
compression and decompression steps followed by energy 
minimization do not target any particular MS packings. 
Yet, the probabilities of different MS packings vary by 
many orders of magnitude. Moreover, as we will show 
below, even when we significantly alter the protocol used 
to generate the MS packings, the most frequent packings 
remain frequent and the rare packings remain rare. 

The question of what gives rise to the extremely diverse 
probabilities of different MS packings is not only impor- 
tant for the analysis of random close packing, but also 
for developing statistical descriptions of granular mate- 



rials and other athermal systems. For example, it has 
been assumed in Edwards' entropy descriptions of gran- 
ular media [33,13 ^h^^ different jammed packings occur 
with approximately equal weights. Similar assumptions 
are also usually employed in PEL theories of glassy ma- 
terials 35]. However, these assumptions should be reex- 
amined, since analogous geometric features are likely to 
control the frequency distributions of relevant states in 
these systems as in our studies of the frequency of MS 
packings. Understanding the reason for the enormous 
probability differences and the comparative roles of the 
frequent and infrequent states are thus important unre- 
solved problems. 

This paper is organized as follows. In Sec. ^1 we de- 
scribe the model and protocol we use to generate the 
MS packings and how we classify and count each distinct 
packing. In Sec. lIIII we review our results for the density 
of MS packings and their frequency distribution in small 
2d bidisperse systems. In Sec. lIVI we investigate the cor- 
relation between local geometric features of configuration 
space near each MS packing and the MS packing frequen- 
cies by adding thermal fluctuations to the packings. We 
look for structural properties that distinguish between 
frequent and rare MS packings in Sec.0 In Sec. I VII we 
introduce a simple phenomenological model to explain 
the dramatic variation in MS packing frequencies. We 
conclude and briefly discuss future research directions in 
Sec.rVIll 



II. GENERATION AND CLASSIFICATION OF 
MECHANICALLY STABLE DISK PACKINGS 

A. Model 

We study the mechanical and statistical properties of 
MS packings in two-dimensional bidisperse systems of N 
frictionless disks interacting via the finite-range, pairwise 
additive, purely repulsive spring potential of the form 



(1) 



Here e is the characteristic energy scale, r.y is the sep- 
aration between particles i and j, Gij — {ai + aj) /2 is 
their average diameter, and Q{x) is the Heaviside step 
function. The particles are in a square unit cell with pe- 
riodic boundary conditions. The mass m of all particles 
is assumed to be equal. 

We consider 50-50 binary mixtures of large and small 
particles with diameter ratio 1.4. We study bidisperse 
mixtures because the presence of particles with different 
sizes inhibits crystallization (provided that the size ratio 
is incommensurate with a binary crystal structure). In 
contrast, monodisperse disk packings crystallize readily 
[36L l3^ , so these systems have entirely different proper- 
ties than glass-forming liquids that are the subject of our 
investigations. Bidisperse mixtures of disks with diame- 
ter ratio 1.4 have been used in several previous studies 
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[Hill m m nil. Such mixtures thus constitute a 
convenient reference system for investigations of funda- 
mental properties of amorphous and glassy materials. 

In systems with finite-range purely repulsive interpar- 
ticle potentials there are two general classes of potential 
energy minima. First, the system can possess a connected 
network of particle overlaps that spans the system; the 
sum of forces on each particle in such networks is zero. 
These configurations have positive total potential energy, 
and each displacement of the particles in the network re- 
sults in an energy increase. For the second type there 
are no particle overlaps, all forces are zero, the total po- 
tential energy vanishes, and the particles can be moved 
without creating an overlap. 

Our focus here is on states that are on the border be- 
tween these two general classes: we assume that the sys- 
tem is at an energy minimum with infinitesimal parti- 
cle overlaps (thus, the total energy is also infinitesimal). 
Since these states are assumed to be in stable mechan- 
ical equilibrium and possess vanishingly small overlaps, 
we term these states MS packings. 

If all particles of a MS packing participate in the sys- 
tem spanning force network, no particle displacement is 
possible without creating an overlap. Occasionally, a 
small number of particles in a MS packing (no more than 
three for the small systems considered herein) do not 
participate in the force network. States with such free 
particles (rattlers) have to be treated with care in the 
packing generation procedures. However, we do not find 
that the properties of MS packings with rattlers deviate 
significantly from those of MS packings without rattlers. 



B. Packing generation protocols 

1. Simulation algorithms 

We use here a class of packing-generation protocols 
that involve successive compression or decompression 
steps followed by energy minimization |l7t l3lj | . The sys- 
tem is decompressed (or, equivalently, the particle diam- 
eters are decreased) when the energy of the system at a 
local minimum is nonzero; otherwise, the system is com- 
pressed. The increment by which the particle packing 
fraction is changed at each compression or decompres- 
sion step is gradually decreased. After a sufficiently large 
number of steps, a MS packing with infinitesimal overlaps 
is thus obtained. 

For each independent trial, we begin the process by 
choosing random initial positions for the particles at 
packing fraction 0o — 0.60 (which is well below the min- 
imum packing fraction at which MS packings occur in 
2d). We then allow the system to relax, and perform 
a sequence of compression/decompression and relaxation 
steps. We can repeat this process for many indepen- 
dent initial conditions, generate a large number of MS 
packings, and measure their respective probabilities for 
a given protocol. 



In the present study, we use two energy-minimization 
methods: (a) conjugate-gradient (CG) minimization al- 
gorithm or (b) molecular dynamics (MD) with dissi- 
pation proportional to local velocity differences. The 
conjugate-gradient method is a numerical scheme that 
begins at a given point in configuration space and moves 
the system to the nearest local potential energy minimum 
without traversing any energy barriers ,3i^]. In contrast, 
molecular dynamics with finite damping is not guaran- 
teed to find the nearest local potential energy minimum 
since kinetic energy is removed from the system at a fi- 
nite rate. The system can thus surmount a sufficiently 
low energy barrier. In the molecular dynamics method, 
each particle i obeys Newton's equations of motion 
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where di is the acceleration of particle i, Vij is the relative 
velocity of particles i and j, fij is the unit vector con- 
necting the centers of these particles, and b is the damp- 
ing coefficient. In the present study, we chose the di- 
mensionless damping coefficient b = ab/y^me = 0.5, but 
this will be varied in subsequent studies. In the infinite- 
dissipation limit 6 — > oo the potential energy cannot in- 
crease during a molecular-dynamics relaxation, and thus 
the molecular-dynamics and conjugate-gradient methods 
should give very similar results. We note, however, that 
even in this limit the two methods are not equivalent 
because there may be more than one energy minimum 
accessible from a given point in configuration space with- 
out traversing an energy barrier. In our previous studies 
[itI I2M I29I I. we used only the conjugate-gradient mini- 
mization algorithm. 



2. Implementation details 

In specific implementations of our MS-packing gener- 
ation protocols, one needs to use appropriate numerical 
criteria for stopping the energy minimization process ei- 
ther in a configuration with overlapping particles form- 
ing a force network or in a state with no particle over- 
laps. For the the conjugate gradient method, we ter- 
minate the minimization process when either of the fol- 
lowing two conditions on the potential energy per par- 
ticle V is satisfied: (a) two successive conjugate gra- 
dient steps n and n -I- 1 yield nearly the same energy 
value, {Vn+i - Vn)/Vn < S = 10-1^; or (b) the poten- 
tial energy per particle at the current step is extremely 
small, Vn < Vnun ~ 10^^® (whcrc V is normalized by 
the energy-scale parameter e). Since the potential en- 
ergy oscillates in time in the molecular dynamics method, 
condition (a) is replaced by the requirement that the rel- 
ative potential-energy fiuctuations satisfy the inequality 

{{v-{v)ry/y{v)<5. 

Following the energy minimization, we determine 
whether the system should be compressed or expanded. 
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If, V < Voiin, the system is below the jamming threshold, 
and thus it is compressed by A(f>. If, on the other hand, 
V > Vmax = 214iin, the system is decompressed by A0. 
For the first compression or decompression step we use 
the packing- fraction increment A<j) — 10"^. Each time 
the procedure switches from expansion to contraction or 
vice versa, Acj) is reduced by a factor of 2. 

For the molecular-dynamics procedure to be efhcient, 
rattler particles with no contacts must be treated with 
care. When the system is near the jamming threshold, 
we set the velocities of rattlers to zero; we also shift the 
center-of-mass velocity of the remaining non-rattler par- 
ticles to zero to assure momentum conservation. This 
modification of the energy-minimization procedure allows 
our stopping criteria to be implemented without change 
even when rattlers are present. Otherwise, the kinetic 
energy of a rattler decays too slowly, and it is extremely 
difficult for the system to reach the threshold Vmin- 

The MS packing generation process terminates when 
Vmin < V < Vmax after energy minimization. In the 
final state, the system is thus in mechanical equilib- 
rium with extremely small overlaps in the range 10^^ < 
1 — Tij/aij < 10~^. We verify the stability of each final 
equilibrium configuration by calculating the dynamical 
or Hessian matrix (i.e., the matrix of second derivatives 
of the total potential energy with respect to the parti- 
cle positions) liO]. In a small percentage of cases we 
find that the dynamical matrix has extra zero eigenval- 
ues that do not correspond to rattlers, which indicates 
that the system is at a saddle point rather than at an 
energy minimum. Such unstable packings are not con- 
sidered. A more detailed discussion of the dynamical 
(Hessian) matrix is given in Sec. Ill CI 

Our procedure for finding mechanically stable disk 
packings allows us to determine the jamming threshold 
in cj) to within 10~^. The procedure is similar to those 
implemented recently for static granular packings with 
and without friction [soIIbH . Our results, however, have 
much greater precision. High accuracy is important in 
our packing-enumeration studies, because some of the 
distinct MS packings have nearly the same packing frac- 
tion. 

To illustrate typical behavior of the system during 
the MS-packing generation process in the version with 
molecular-dynamics energy minimization procedure we 
show, in Fig. |21 the evolution of the potential energy per 
particle V for several consecutive values of (j). Note that 
the potential energy is not monotonic in time because of 
finite particle inertia. 



C. Identification of Distinct MS Paclcings 

We consider two MS packings to be identical if they 
have the same network of contacts (i.e., their networks 
can be mapped onto each other by translation, rotation, 
or by permutation of particles of the same size). An 
analysis of the topology of the network of contacts is. 
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FIG. 2: Potential energy per particle V as a function of time 
t at several (p during the process of creating a mechanically 
stable packing of 10 particles at (j>c ~ 0.8244438. The curves 
correspond to (j) ^ 0.8245 (dotted), 0.82443 (dashed), 0.82444 
(long-dashed), 0.824443 (dot-dashed), and 0.8230 (solid). MD 
energy minimization was employed. Note that for (j> > 4'c 
energy minimization terminates a.t V > Vmin = 10~^^, while 
for (j> < (j)c minimization terminates when V — Vmin- The 
spring timescale ts = a^frnfk, where a is the small particle 
diameter, was chosen as the normalization for time. 

however, a fairly complex task. Thus, in practice it is 
convenient to use some alternative but equivalent crite- 
ria. 

A test based on the packing fraction alone is inade- 
quate, because some packings with distinct topological 
networks have the same packing fraction (within the nu- 
merical noise). We find, however, that a test based on 
a comparison of the spectra of the dynamical matrix is 
sufficient. 

For a pairwise additive, rotationally invariant poten- 
tial Q the dynamical matrix (Hessian) is given by the 
expressions 5] 



M; 



"(^Q/3 ^ija^ijl3^ CijTijct'l^ijf^ J i 7^ J, (3) 



and 



E 



M, 



(4) 



where — dV/drij and 



d'^V/dr 



In the above 

relations the indices i and j refer to the particles, and 
a, 13 = x,y represent the Cartesian coordinates. For a 
system with Nf rattlers and N' ~ N — Nf particles form- 
ing a connected network the indices i and j range from 
1 to N', because the rattlers do not contribute to the 
potential energy. 

The dynamical matrix is symmetric, and it has dN' 
rows and columns, where d = 2 is the spatial dimen- 
sion. Thus it has dN' real eigenvalues {m^}, d of which 
are zero due to translational invariance of the system. 
In a mechanically stable disk packing, no set of particle 
displacements is possible without creating an overlapping 
configuration; therefore the dynamical matrix has exactly 



6 



TABLE L Number of distinct MS packings nf° and 
obtained and number of trials and n^'^ performed using 
the MD and CG energy minimization methods and estimated 
total number of distinct MS packings nl°^ calculated using an 
extrapolation of the CG results for several system sizes A'^. 
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d{N' — 1) nonzero eigenvalues. In our simulations we use 
the criterion \mi\ > m^nin for nonzero eigenvalues, where 
Wniin = 10~^ is the noise threshold for our eigenvalue 
calculations. 

In our numerical simulations we distinguish distinct 
mechanically stable disk packings by the lists of eigen- 
values of their dynamical matrices. We assume that two 
MS packings are the same if and only if they have the 
same list of eigenvalues. (The eigenvalues are considered 
to be equal if they differ by less than the noise threshold 
Wniin-) To verify our method we have also compared the 
topology of the network of particle contacts in different 
packings, and we have found that these two methods of 
identifying distinct mechanically stable packings agree. 

As noted above, it is in general not true that each dis- 
tinct MS packing possesses a unique packing fraction 0. 
However, we find that for small systems only at most 
a few percent of distinct MS packings share the same 
packing fraction. Thus, in the following we will asso- 
ciate a unique (fi with each MS packing to simplify the 
discussion. Also, approximately 10% of the distinct MS 
packings contain at least one rattler particle. In these 
configurations, we ignore the translational degeneracy of 
the rattlers — two configurations that have the same con- 
tact networks of non-rattler particles are treated as the 
same. We do not find that the properties of the MS pack- 
ings with rattlers deviate significantly from those of the 
MS packings without rattlers. 



have been employed to determine the dependence of the 
results on the packing-generation protocol. 

In Table IHwc list the numbers of distinct MS packings 
and nj*^ obtained using the MD and CG methods 
and the corresponding numbers of trials performed nf^^ 
and n^*^ . Our estimate for the the total number of MS 
packings that exist n'°* for each system size is also given. 

The total number of distinct MS packings has been 
estimated by extrapolating the relation between n^*^ and 
nf^ using the approach proposed in |l3|. For iV < 10 
the number of distinct packings n^*^ generated by the 
CG method for the given number of trials saturates, 
which indicates that nearly all MS packings have been 
obtained for these system sizes For A*" = 12 and 14 
the CG method has yielded about 90 % and 70 % of the 
total MS packings, respectively. 

For the two smallest systems studied, A^ = 4 and 6, the 
sets of MS packings generated by the CG and MD meth- 
ods are identical. For larger system sizes, the CG method 
finds more packings than the MD method for a fixed 
number of trials because a large fraction of MS pack- 
ings become extremely rare when using the MD method. 
(This interesting feature of the MD packing-generation 
algorithm will be discussed in more detail below.) 

The results in Table |l] indicate that the number of dis- 
tinct MS packings grows exponentially with increasing 
system size. In addition, the number of trials needed to 
find all MS packings (or to find a large fraction of them) 
exceeds by orders of magnitude the number of distinct 
MS packings themselves. For example, with as many as 
3 X 10^ MD trials, we have found only about 1250 pack- 
ings out of approximately 1600. The CG results show 
similar behavior, but rare packings are not as rare as 
for MD. These results indicate that the frequencies with 
which MS packings occur can vary by many orders of 
magnitude. In this article we seek to understand the 
source and the significance of this property. 



B. Frequency distribution and density of MS 
packings 

1. Protocol- dependent and protocol-independent quantities 



III. PROBABILITY DISTRIBUTION OF MS 
PACKINGS 

A. Total number of distinct MS packings 

We have applied our algorithm for finding mechanically 
stable disk packings using a large number of independent 
trials with different starting configurations for systems 
with up to A^ = 256 particles. Here we focus, however, 
only on small systems with A^ < 14, because for these 
values of A^ we were able to generate a significant fraction 
of all distinct MS packings. Both the conjugate gradient 
and molecular dynamics energy minimization techniques 



Nearly complete enumeration of the MS packings for 
small systems allows us to characterize the distribution 
of MS packings in much more detail than is possible for a 
system with large A^. For sufficiently small systems, each 
distinct MS packing can be generated multiple times; 
thus for each distinct packing (indexed by k) , we can de- 
termine the protocol-dependent probability fk for which 
it occurs. Next, for a given packing-fraction interval 
we can evaluate the probability 

Pi^)Aq} = np{Aq})/nt (5) 

of generating a state in the interval A0 for a specific 
protocol. We can also define the protocol-independent 
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FIG. 3: The probability density P(<^) (solid), density of MS 
packings p{(j)) (dotted), and the frequency distribution /3(0) 
(long-dashed) averaged over bins with width A(j} — 0.02 for 
(a) 10- and (b) 14-particle systems using the CG energy min- 
imization method. Note that the results for N = 14 are ap- 
proximate since we have not found all of the MS packings for 
this system size. 

number of distinct MS packings 



p{4>)A(j) = ns(A(/))/ns 



(6) 



that exist in this interval. In Eq. [Sj np{A<p) denotes 
the number of trials that produce MS packings in the 
interval of interest and nt is the total number of trials. 
The corresponding quantities in Eq. are the number 
of distinct MS packings ns(A0) that exist in the given 
packing-fraction interval and the total number of distinct 
MS packings Ug in the system. We note that both the 
probability density P{(p) and the density of MS packings 
p{4)) are normalized to unity. 

In addition to the quantities defined by Eqs. El and El 
we also consider the average frequency distribution 



(7) 



for MS packings in the interval A0 near the packing frac- 
tion (f>. The continuous frequency distribution and the 
discrete probabilities fk for obtaining the fcth packing 
satisiy the relation 



where 



{fk)= E h/MAq^) 



(8) 
(9) 
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FIG. 4: A comparison of the total probability distribution 
P{(j>) (solid) and frequency distribution P{4>) (long-dashed) 
averaged over bins with width A<j) = 0.02 for a 10-particle 
system using the CG (thick lines) and MD (thin lines) meth- 
ods. The density of MS packings (dotted) is also shown. The 
MD results only include approximately 75% of the total MS 
packings. The inset shows /3('/>) using the CG (squares) and 
MD (circles) methods for the same system on a logarithmic 
scale. The slopes of the lines are approximately 25 and 33. 

is the average value of the probability fk in the interval 
A0. 



2. Relation between density p{(j)) of MS packings and 
random-close packing 

The decomposition of the probability density jSl into 
the density of MS packings O and the average frequency 
13 was studied previously in our recent article [131 for 
the CG packing-generation protocol. Sample results of 
our simulations are shown in Fig. j^l where we plot the 
functions P(0), p{(f>), and /3((/)) for = 10 and 14. 

Two key features should be noted from this figure. 
First, the peaks of the density of MS packings p{(t)) and 
probability distribution P{(j>) become narrower as the 
system size grows. Second, the separation between the 
peaks decreases with increasing N. We observe that the 
distance between the maxima of /3((/)) and P((/)) is of the 
order of the peak width both for = 10 and 14. 

This behavior suggests that the random-close packed 
density 0rcp can be re-defined as the position of the peak 
in the density of MS packings as we proposed in Ref. 
[l7| . The fluctuations in the packing fraction of large 
MS packings can be described by a superposition of local 
fluctuations in a large number of subsystems. Therefore, 
the width of the peak in p((f>) should scale approximately 
as 7V~^/^ by the central-limit theorem, which was con- 
firmed in Ref. j^^. The position of the peak in p(0) 
will coincide with the peak position of the probability 
distribution P{(j)) = f3{(j))p{(j)) in the large system limit, 
unless the frequency (3{4>) varies so rapidly with (j) that 
it elevates the tail of the distribution p{4>). 

According to this definition, random close packing is a 
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FIG. 5: The discrete probability to obtain a MS packing at 
a given packing fraction (j>k for N = 10 particles using the MD 
energy minimization method. The inset shows a magnified 
view of the probability on a logarithmic scale over a narrow 
region of packing fraction between 0.81 and 0.82. Even within 
this narrow region, the probability varies by more than five 
orders of magnitude. 

protocol-independent quantity. Our picture is thus con- 
sistent with the observation that for a large class of pro- 
tocols essentially the same random close packed density 
is obtained for macroscopic systems. This picture is also 
supported by the results shown in Fig. 01 where the prob- 
abihty density P(0) is depicted for the CG and MD pro- 
tocols. The results indicate that the peaks in P{(j)) nearly 
coincide for the two protocols, in spite of a noticeable 
difference between the corresponding frequency distribu- 
tions for the MD and CG methods, especially at the 
lower end of the range in packing fraction. 

The above-described scenario seems quite plausible, 
but there are still a number of puzzling findings that 
need to be addressed before amorphous jammed pack- 
ings are fully understood. First, we find that even for the 
class of algorithms considered here the frequency (3{(j)) is 
a rapidly varying, exponential function of 0. Such ex- 
ponential behavior of /?(</)) is clearly seen in the inset of 
Fig. 2] for the 10-particle system. Similar results were 
also obtained for larger N. 

The results shown in Fig. 0] indicate that the exponen- 
tial variation of the frequency distribution f3{4)) is quite 
rapid — we find that this quantity changes by several or- 
der of magnitude in the relatively narrow range of pack- 
ing fractions where a significant number of MS packings 
exists. Moreover, the exponent is an increasing function 
of N, as revealed by simulations presented in Ref. [T^ . 
Under the assumption that the exponential behavior also 
holds for large systems, the peak of the probability P{(j>) 
is controlled by the peak of the density p{(j)) of MS pack- 
ings, provided that \og(5/N — > for iV ^ oo. With 
current computational resources we are not yet able to 
assess the validity of this condition. The source of the 
rapid variation of /3(0) and the system-size dependence 
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FIG. 6: (a) The discrete probabilities /fc for MS pack- 
ings obtained using the CG method within a narrow interval 
A(/!> near the peak in p{(j)) normalized by the maximal prob- 
ability in the interval J™^", sorted in increasing order, and 
plotted on a logarithmic scale for A'^ = 10 (dotted), 12 (long- 
dashed), and 14 (solid) particles. The interval = 0.02 for 
TV = 10 and 12 and 0.004 for iV = 14. The index k on the 
horizontal axis is normalized by the total number of distinct 
MS packings Us that exist within A(/>. (b) The probabili- 
ties logj^Q ifk/fk^^'^) scaled by an 0(1) constant a. a and Us 
were chosen to yield the best collapse of the scaled curves for 
A'^ = 12 and 14 with the unsealed curve for A*' = 10. 



of this variation are thus important unresolved issues. 

Our simulations also reveal another striking result. We 
find that the discrete probabilities of MS packings can 
differ by many orders of magnitude even when they pos- 
sess nearly the same value of (j). This behavior will now 
be discussed. 



3. Discrete probabilities of distinct MS packings 

The discrete probabilities fk of distinct MS packings 
in the range of packing fractions within the peak of the 
probability density P{(j)) are depicted in Fig. for the 
10-particle system. The results are shown for the MD 
packing-generation protocol; the corresponding distribu- 
tion for the CG protocol is similar. 

Consistent with our findings for the continuous fre- 
quency distribution (3{(t>) the probabilities fk are, on av- 
erage, significantly larger at larger values of (f). The dis- 
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FIG. 7: (a) A comparison of the sorted discrete probabilities 
fk normalized by J^'^'' for the CG (soHd line) and MD (dot- 
ted line) methods for an interval in packing fraction near 
the peak in p{(j)) for A*' = 10. (b) The discrete probabilities 
fk for the 20 most frequent MS packings depicted in (a) from 
the CG packing-generation method (filled squares) are com- 
pared to the probabilities for the same packings obtained from 
the MD method (fiUed circles). A similar comparison of the 
probabilities for 25 less frequent packings obtained from the 
CG (open squares) and MD (open circles) methods are also 
shown. The circled regions in (a) identify the MS packings 
that were compared in (b). 



Crete probabilities /fc, however, also vary dramatically 
from one MS packing to another. They can differ by 
more than five orders of magnitude even in a narrow in- 
terval of (j), as illustrated in the inset of Fig. |S1 These 
large local probability fluctuations occur over the entire 
range of (f>. 

It is clear from these results that MS packings do not 
occur with approximately equal probabilities as has been 
assumed in the Edwards' entropy descriptions of powders 
and granular media 33, 34] . The large variation of the 
probabilities is a rather puzzling result because our algo- 
rithms do not target specific packings. This suggests that 
the probabilities are determined by geometrical features 
of the energy landscape. We will return to this important 
problem in Sec. IIVI but let us first examine the discrete 
probabilities fk in more detail. 

Figure shows the discrete probabilities fk for MS 
packings within a single packing-fraction interval Acf) 
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FIG. 8: (a) A comparison of the sorted discrete probabili- 
ties fk normalized by /max for the CG (solid line) and MD 
(dotted line) methods including all MS packings for N = 10. 
/max is the frequency of the most probable MS packing for 
A'' = 10. (b) The discrete probabilities fk for the most fre- 
quent MS packings depicted in (a) from the CG packing- 
generation method (filled squares) are compared to the prob- 
abilities for the same packings obtained from the MD method 
(filled circles). A similar comparison of the probabilities for 
less frequent packings in the packing-fraction range [0.76, 0.82] 
obtained from the CG (open squares) and MD (open circles) 
methods are also shown. The circled regions in (a) identify 
the MS packings that were compared in (b). 

near the maximum of the density of MS packings p(<^) 
for the CG method. The probabilities fk are sorted in 
increasing order and are plotted versus the index k nor- 
malized by the total number of states ^^(At/)) that ex- 
ist in the interval. The probabilities are normalized by 
the maximal probability value f^'^^ (within A0) and dis- 
played on a logarithmic scale. Results for systems with 
N — 10, 12, and 14 particles are provided. 

It is important to notice that all three curves in Fig.|^ 
have a similar shape. Indeed, by rescaling the data 



iogio(/fc//r'^)-«iogio(/fc//r'^) 



(10) 



by a factor a — 0{\) the results for different system sizes 
collapse onto a single master curve, as depicted in Fig. 
El (b) . Since the total number of MS packings in these 
three systems differs by more than two orders of magni- 
tude, the scaling^] is a non-trivial result. This result 
further reveals that there must be geometric features of 
configuration space that determine the packing probabil- 
ities. 
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FIG. 9: (a) The truncated density of only the most frequent 
MS packings Pq{(j)) for the CG energy minimization method 
in a 10-particle system. The accumulated probability in each 
packing fraction interval A(j) = 0.02 is given by the param- 
eter q; we display q — 0.3 (dotted), 0.5 (dashed), 0.7 (long- 
dashed), 0.9 (dot-dashed), and 1.0 (solid), (b) The fraction of 
distinct MS packings nq{<j)) that contribute to Pq{4>) in each 
<j!> interval for the same system in (a). 

We note that for N = 12 and 14 the total number of 
distinct packings in the interval A(f) has not been mea- 
sured directly — the most infrequent states could not be 
generated because we were not able to perform a suffi- 
ciently large number of trials with our current computa- 
tional resources. The estimates for ns{A(j)) used to scale 
the horizontal axis in Fig.|ni(and subsequent figures) were 
obtained by matching the rescalcd discrete probabilities 
for = 12 and 14 to the corresponding curve for = 10, 
for which the total number of MS packings is known. 

Unlike the density of MS packings, the packing proba- 
bilities fk are protocol dependent. To examine the effect 
of different protocols on the probability distribution we 
compare, in Fig. [71 a), the sorted probabilities fk (scaled 
by f^^^) for the CG and MD protocols in a 10-particle 
system over a small interval in packing fraction near 
the peak in p{4)). We observe that the relative proba- 
bilities are significantly lower for the MD protocol, ex- 
cept in the high-probability regime, where both sorted- 
probability curves coincide. 

The sensitivity of individual packing probabilities to 
the change of the packing-generation protocol is exam- 
ined in Fig. EI&). In this figure we compare the CG and 
MD probabilities of the 20 packings that occur most fre- 



quently in the interval Acf) according to the CG proto- 
col. We also show 25 packings that have much lower 
frequencies. These results indicate that while the indi- 
vidual probabilities may change even by several orders of 
magnitude (especially in the low-probability regime), a 
significant shuffling between the sets of frequent and in- 
frequent packings does not occur when the energy min- 
imization protocol is changed. The packings that are 
frequent according to the CG protocol typically become 
even more frequent according to the MD method, and 
the states that are infrequent become even less frequent. 
Moreover, the sets of the most frequent packings for the 
two packing-generation methods nearly coincide. 

Similar observations also hold for the set of all MS 
packings (rather than only those within a small interval 
A0), as can be seen from the results shown in Fig. |S1 In 
particular, we find that for the 10-particle system ^ 80 
out of the 100 most frequent CG packings are also most 
frequent when they are generated using the MD protocol. 

We believe that the above results profoundly affect the 
way one should interpret and explain a range of jamming 
and glassy phenomena including random-close packing, 
dense slow granular flows, and arrested dynamics in glass- 
forming liquids. For example, the Edwards' entropy de- 
scription of nearly jammed granular materials is based 
on the assumption that different jammed packings occur 
with approximately equal weight. Our findings indicate, 
however, that this assumption should be reconsidered. 

Should we include all of the MS packings in statis- 
tical descriptions of jammed and glassy systems even 
though the probabilities vary so strongly or should we 
include only the most frequent packings? If only the 
most frequent packings are needed to describe jamming 
and glassy phenomena, how do we distinguish between 
the frequent and infrequent ones? Similar questions ap- 
ply to not only to Edwards' entropy descriptions, but also 
to definitions of random close packing and descriptions of 
glassy materials in terms of local minima in the potential 
energy landscape. 

Results from preliminary investigations aimed at ad- 
dressing these questions are shown in Fig. |51 In panel 
(a), we plot the truncated density of MS packings Pq{<j)), 
which is defined as the density of only the most frequent 
MS packings in each interval of 0. The accumulated 
probability of these packings (normalized to unity in each 
packing-fraction interval) is given by the parameter q. 
For future use, the set of all states contributing to pq is 
denoted by {Vk}'^- The results in Fig. |5| are displayed 
for the CG packing-generation protocol and iV = 10. 

Consistent with the simulation data presented in Fig. 
I^lfor an interval near the maximum of the density of MS 
packings, the most frequent 15% of MS packings con- 
tribute as much as 50 % of the local probability, accord- 
ing to the results in Fig. 0(b). Also, to accumulate 90% 
of the probability we only need roughly 50 % of the MS 
packings. 

One of our key observations is that the density of MS 
packings is insensitive to the cutoff probability q. Thus, 
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the peaks of Pq{(j)) for different probability-truncation 
levels nearly coincide. In contrast, the peaks of 
and P{(j>) in Fig. E^a) are more separated (although still 
within the peak width), because the average frequency 
distribution depends exponentially on (j). 

The simulation results described so far reveal that the 
density of discrete MS packings is an important quantity 
that may control, for example, the position of the peak in 
the probability distribution of MS packings for large sys- 
tems, independent of the compaction protocol. However, 
we have also found that the probabilities of individual 
MS packings vary dramatically from one packing to an- 
other. We have also discovered interesting regularities of 
the packing distributions, both as a function of 4> a-nd 
within narrow packing-fraction intervals. 

We believe that our observations are important for con- 
structing theoretical descriptions of amorphous packings 
and slow dynamics in glassy materials. Perhaps only the 
most frequent MS packings control the behavior in these 
systems. However, we must first understand what deter- 
mines the frequency with which MS packings occur in or- 
der to identify correctly the set of frequent MS packings. 
In the following sections we investigate which structural 
and geometric factors play an important role in deter- 
mining the MS packing frequencies. 



IV. POTENTIAL ENERGY LANDSCAPE NEAR 
MS PACKINGS 

In this section we examine whether or not the proba- 
bilities of MS packings can be correlated with local char- 
acteristics of the potential energy landscape in the neigh- 
borhood of each packing. The local features include the 
heights of the energy barriers separating different MS 
packings and the shape and size of the local region of con- 
figuration space that is visited when mechanically stable 
packings are subjected to thermal fluctuations. As de- 
scribed below, we find interesting correlations between 
the probabilities of the MS packings and these local geo- 
metric features. Our results, however, show that a more 
global analysis of the topography of the potential energy 
landscape is also required. 

As in our packing-generation procedures, we consider 
a system interacting via the one-sided soft spring poten- 
tial^ However, our results can be applied more broadly. 
First of all, local fluctuations around a MS packing in 
a system interacting via a finite-range repulsive poten- 
tial can approximately be mapped onto the motion of 
a hard-disk system. The effective disk radii correspond 
to the distance at which the pair potential V{rij) equals 
the thermal energy T Thus, adding thermal fluctu- 
ations to a MS packing is analogous to decompressing a 
collectively jammed hard-disk system and then perform- 
ing hard-particle energy-conserving collision dynamics. 

Second, crossing energy barriers and transitioning from 
one MS packing to another resembles the evolution from 
the basin of one inherent structure 1431 to another in 
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FIG. 10: The packing fraction (j> of the MS packing that 
is nearest to the instantaneous MD configuration at time t 
(normalized by tkin) for temperatures T = 2.5 x 10~® (solid 
line) and a factor of two lower (dotted line). The minimal 
breaking times r for these two temperatures are indicated by 
the arrows. 

glass-forming systems that interact via continuous poten- 
tials. We note that there are many more inherent struc- 
tures than MS packings at each so that a complete enu- 
meration becomes prohibitively costly even for extremely 
small systems, unlike enumeration of the MS packings. 
Also, MS packings can be interpreted as metabasins |4J| 
of the PEL for systems with finite-range repulsive poten- 
tials. Thus our results are important not only for under- 
standing random-close packing but also glassy behavior 
in soft- and hard-particle systems. 



A. Breaking times and energy barriers 

1. Measurement of breaking times 

In our approach, we probe the local potential energy 
landscape near a given MS packing via thermal fluctua- 
tions. Each initial MS-packing configuration is thermally 
excited by adding kinetic energy to the system. The en- 
ergy is introduced by choosing the initial velocities of 
the particles randomly from a Gaussian distribution with 
variance 2T. We then allow the system to evolve at con- 
stant energy according to the evolution equation [5] with 
the damping coefficient b set to zero. Even though the 
systems we study are small, they behave like thermal 
systems in the sense that the 2T energy input is quickly 
partitioned equally among the configurational and kinetic 
degrees of freedom: {V) ^ T and {K) ~ T, where K is 
the kinetic energy per particle. 

During the course of the molecular-dynamics run, we 
periodically save the particle coordinates. For each 
snapshot, the particle coordinates are fed into our MS 
packing-generation routine to find the nearest MS pack- 
ing using the CG energy minimization scheme. The near- 
est mechanically stable packings for each snapshot are 
then compared to the original MS packing at t = 0. 
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This procedure allows us to measure the first-passage 
breakup time r for the system to make a transition from 
the original MS packing to a different one. In this section, 
all times are measured in units of the kinetic timescale 
ikin = cr-y/ m/T. We chose the time interval between 
instantaneous MD snapshots to be 5-20 At ^ r (where 
At is the integration time step) , which is small enough so 
that it does not influence our results. For each initial MS 
packing and temperature T, the procedure is repeated at 
least 20 times with random initial velocities to obtain the 
average breaking time (r) . 

A typical example of the breaking-time measurements 
at two different temperatures for a MS packing near the 
peak of the density of states p[<f)) is depicted in Fig. ^1 
The figure shows the evolution of the packing fraction (p 
of the MS packing that is nearest to each instantaneous 
MD configuration. After the state breaks away from 
the neighborhood of the initial MS packing, it cascades 
through a set of MS packings with gradually increasing 
4> and ends up oscillating between several high-packing- 
fraction, high-probability MS packings. We emphasize 
that the molecular dynamics evolution takes place at con- 
stant <f>. The separate compression/decompression and 
energy-minimization procedure is used only to identify 
the corresponding MS packing for each instantaneous 
configuration from the MD evolution. 

Our method for finding the nearest MS packings is sim- 
ilar to thermally quenching instantaneous MD configura- 
tions to the nearest inherent structure or local potential 
energy minimum However, apart from energy mini- 
mization, in our procedure we incorporate the additional 
steps of shrinking and growing the particles to arrive at 
a MS packing with infinitesimal particle overlaps |45| . 

2. Measurement of energy harrier tieights 
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FIG. 11: The average breaking time (r) required for the 
system to surmount the lowest energy barrier eo versus inverse 
temperature 1/T for three different MS packings with A'^ = 10 
particles. The slope of the solid lines in panels (a)-(c) are 
eo PS 3.3 X 10"^ 1.4 x 10"^ and 1.1 x 10"". 



In order to determine the height of the energy barrier 
that separates a given MS packing from other packings, 
we have performed a series of breaking-time measure- 
ments over a wide range of temperatures T . The simu- 
lations were performed on a randomly selected group of 
105 MS packings with N = 10. 

Measurements of the average minimal breaking time 
(r) versus the inverse temperature for three different MS 
packings are presented in Fig. Each data point in 
these plots was obtained from an average of at least 20 in- 
dependent measurements of the breaking time; the stan- 
dard deviation of r is comparable to the symbol size. 

The results in Fig. II II indicate that for sufficiently low 
temperatures the system exhibits Arrhenius behavior 

(r) = Tooe^"/^. (11) 

The quantity eg in the above relation is interpreted as 
the height of the energy barrier that separates a given 
MS packing from others. Below, both eq and T will be 
measured in units of the characteristic energy scale e of 
the repulsive spring potential in Eq.^ The parameter Too 



characterizes the timescale at which the system explores 
the energy landscape. 

We find from our results that the exponential behav- 
ior is quite robust despite the fact that we use small 
systems with only ^ 20 translational degrees of freedom. 
The evolution occurs at constant energy with no heat 
reservoir — yet we observe pronounced thermal behavior. 
We note that the breaking time in the Arrhenius regime 
varies by several orders of magnitude, which allows us to 
determine the barrier height eo quite accurately. 

At high temperatures, the MS packings break to a 
number of different destination MS packings because the 
system possesses enough kinetic energy to traverse a 
broad range of barriers. In the low-temperature regime, 
the system breaks by jumping over the lowest energy bar- 
rier, unless there are several barriers of nearly the same 
magnitude. By fitting the simulation data (such as those 
presented in Fig. Illfl to the Arrhenius form ^] in the 
low-temperature regime we can thus measure the low- 
est energy barrier eo for each MS packing studied. At 
moderate temperatures, we can separately calculate the 
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FIG. 12: The discrete probability fk with which a given MS 
packing occurs versus the lowest energy barrier eo associated 
with that packing using the (a) CG and (b) MD energy mini- 
mization methods. All 105 MS packings in our sample are in- 
cluded. The filled circles correspond to the most probable set 



of MS packings {Vk} 



while the open circles indicate the 



remaining less frequent ones. The long-dashed line in panel 
(b) with slope 1.7 points out the correlation between and 
eo. The open square (cross) in both panels corresponds to the 
twin packing (meta-packing) discussed in Sec. IIV A 31 

average breaking time (r) for each destination MS pack- 
ing, which, in principle, allows us to measure the lowest, 
next lowest, and subsequent higher energy barriers. 



3. Relation between energy barriers an probabilities of MS 
packings 

The Arrhenius plots presented in Fig. reveal that 
the magnitudes of the minimal energy barriers cq sepa- 
rating distinct MS packings vary by many orders of mag- 
nitude. It is thus interesting to determine whether or not 
there is a relation between the energy-barrier heights and 
the MS packing probabilities /fc, which exhibit a similar 
large variation, as discussed in Sec. Illll We have mea- 
sured the energy-barrier heights for a sample of 105 MS 
packings, which comprise only about 6.5% of the total 
number of MS packings for the 10-particle system con- 
sidered (cf. the results in Table . The small sample 
size does not allow for a complete quantitative analysis 
of the problem, but we are able to make some interesting 
qualitative observations. 

In Fig. 1121 we compare the minimal energy barrier 



heights Co and the corresponding probabilities fk for our 
sample of MS packings. The probabilities shown in panel 
(a) were obtained using the CG protocol, and those de- 
picted in panel (&) correspond to the MD algorithm. In 
both plots, the filled symbols indicate the data points 
for the packings in the set {Vk}'^ of the most probable 
packings contributing to the truncated density of states 
Pqicj)) at the probability truncation level q = 0.70. (The 
set {Vk}'^ is defined near the end of Sec. 11116 31 1 The 
open symbols correspond to less probable packings. 

There are several interesting features to be noticed in 
these plots. In the case of the CG protocol, illustrated in 
Fig. ll2r a'l. we see only a gross correlation between fk and 
eq. The high- probability states tend to have high energy 
barriers, and the barriers of low-probability packings do 
not exceed 10~^ in our sample. Otherwise, the scatter in 
the data is very large; for example, the energy barriers of 
the packings in the probability range 10~^ to 10"'^ vary 
by nearly ten orders of magnitude. 

In contrast, there is a much stronger correlation be- 
tween the probabilities and energy barriers in the subset 
of packings that have the maximal probability in each 
local packing-fraction region. We note that for a given 
value of fk, the packings from the high-probability sub- 
set {T-'fc}''^^ '' (shown as filled circles in Fis. \V2i a) ) tend 
to have the highest energy barriers eo for a given fk- 

The results presented in Fig. \l2i b) show that the 
points in the high-probability subset are insensitive to 
the packing-generation protocol, consistent with the re- 
sults shown in Fig.^ However, the packings outside the 
subset {'Pfc}''^°'^ shift significantly in the low-probability 
direction when switching from the CG to MD method. 
Since packings outside {Pfe}'^^^ '' tend to have low energy 
barriers, the MD probabilities are much more strongly 
correlated with the heights of the minimal energy bar- 
riers than the CG probabilities. In fact, the data for 
the most probable states {7-"^}'^°'^ in Fig. ll2r &') roughly 
scale as a power-law ^ Eq with A ~ 2. 

The above observations point out that the relationship 
between the probabilities of MS packings and eo is com- 
plex. On the one hand, the magnitude of the minimal 
energy barrier separating a given MS packing from other 
packings is not directly linked to the packing probability 
fk, especially for the CG protocol. On the other hand, 
the energy-barrier heights clearly determine some impor- 
tant features of the probability distribution. In particu- 
lar, the probabilities fk of MS packings with low energy 
barriers depend strongly on the packing-generation pro- 
tocol, which can be seen by comparing the results in Figs. 
El (a) and (&). Moreover, there is a significant correla- 
tion between eo and fk for a subset of the most frequent 
packings for each (f). 

Strong protocol dependence of the probabilities fk for 
MS packings with low values of eo can be understood 
intuitively by noting that a system with weak but non- 
vanishing thermal fluctuations is able to surmount low 
energy barriers but not high ones. Thus for protocols 
such as the MD algorithm, the low-barrier packings tend 
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FIG. 13: The packing fraction <j) of the MS packing that is 
nearest to the instantaneous MD configuration at time t for 
(a) a temperature comparable to the minimal energy-barrier 
height T ~ eo and (b) a temperature that is four times larger. 

to have low probabilities. After traversing the low barri- 
ers, the system becomes trapped in packings with higher 
barriers, which thus have higher probabilities. 

We expect similar behavior to occur during the relax- 
ation process of glass-forming liquids following a thermal 
quench. During the quench, the system first quickly re- 
laxes by surmounting low energy barriers. At later times 
it evolves much more slowly through a set of high-barrier 
states, which is closely related to the high-probability 
subset {"Pfc}^ of MS packings. During the slow-evolution 
stage, the properties of the system are thus strongly af- 
fected by the truncated density of the most frequent 
states similar to the truncated density of MS packings 
Pq{4>) introduced in Fig. |3 

The large scatter of the data points in Fig. ^] indi- 
cates that additional geometrical features of the poten- 
tial energy landscape — not simply the minimal energy 
barriers — must influence the probabilities of MS pack- 
ings. Below we qualitatively discuss a few nonlocal topo- 
graphic features that may play a significant role. A more 
detailed, quantitative investigation of such features is be- 
yond the scope of the current study, but will be pursued 
in future work. 

Twin packings and meta-packings A relatively sim- 
ple feature of the PEL that can give rise to significant 
fluctuations in the relation between the minimal energy 
barrier eg and the probabilities fk is illustrated in Fig. 



[T51 This figure shows two sample breakup trajectories 
for the MS packing indicated by the open square in Fig. 
rr?] The trajectory depicted in Fig. [T^ a) corresponds to 
an initial thermal excitation with T comparable to the 
minimal barrier height eo; the trajectory depicted in Fig. 
I13r 6) corresponds to a thermal excitation that is four 
times larger. 

According to the results in Fig. E| at the lower tem- 
perature the system oscillates between two MS packings 
that have nearly the same packing fraction cf). We refer to 
these as twin packings. In order for the system to make 
a transition to a non-twin packing within this timescale, 
the temperature must be increased dramatically. After 
surmounting the barrier that takes the system to a non- 
twin packing, the system undergoes the usual cascade of 
transitions through a sequence of packings with increas- 
ing <j) (cf. the trajectory depicted in Fig. I10|l . 

The meta-packing that consists of the two twin pack- 
ings between which the system oscillates at low temper- 
atures can have a much higher energy barrier than the 
barrier cq separating the two component packings. We 
expect that the probability of the twin packings is con- 
trolled by the meta-packing energy barrier ei rather than 
the minimal barrier eo <C ei. Thus, we have also included 
the data point (ei, fk) for the twin-packing in Fig. 1121 us- 
ing a cross as the symbol. Note that using ei instead of eo 
brings this point closer to the cluster of data for the most 
probable subset of packings {Vk}'^^'^ '^. A generalization 
of this picture to a larger group of p ackings will be car- 
ried out in a future publication |4y|. When we perform 
this analysis we may likely find that more than two com- 
ponent packings can belong to a given meta-packing. We 
note that the concept of meta-packings is closely related 
to that of metabasins of the PEL P, |3| . 

Probability streams Identifying neighboring packings 
that are separated by small energy barriers and grouping 
them into meta-packings may decrease the scatter of the 
data in Fig.^l but it is unlikely to eliminate it. Thus, it 
is apparent that not only energy barriers but other geo- 
metric features of the PEL are important for determining 
the probabilities fk- At present, we do not have any di- 
rect measurements of these additional features. However, 
our simulation results provide important clues that indi- 
cate directions for further investigations of the problem. 

We may consider here two possibilities. First, one 
could assume that the volume Qk of the local region in 
configuration space with energy y < eo in the neighbor- 
hood of a given MS packing determines its probability 
fk- For a given eo, this volume may have large varia- 
tions due to variations of the shape of the local potential- 
energy basin from one MS packing to another. However, 
as shown in Sec. IIV 51 below, this possibility can be ex- 
cluded since the multi-dimensional volume flk depends 
too strongly on the height of the energy barrier. 

Thus, a more complex, non-local scenario may better 
explain the MS packing probabilities. According to this 
scenario, the probability fk of a MS packing depends on 
two competing factors: the probability of getting into 
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the neighborhood of that packing and the probabihty of 
leaving the local region before the MS packing is reached 
during the relaxation process. The probability of leaving 
the local region is largely controlled by the height of the 
energy barrier, and thus it is strongly protocol depen- 
dent. The probability of arriving into the neighborhood 
of an MS packing is affected by the size of the local region 
(therefore states from the high probability subset {Vk}'^ 
tend to have high energy barriers) , but it is more strongly 
influenced by the chain of events that brought the system 
into the neighborhood of the MS packing. These events, 
in turn, are determined by the features of PEL that can 
be far away from the local potential-energy basin for a 
given MS packing. 

The large scatter in the probabilities for a given value 
of the energy barrier indicates that the flux of probability 
during the compaction process is very nonuniform. It ap- 
pears that there are many "dry" regions, with very small 
probability flux, and "probability streams," where the 
probability flux is very large. If a given MS packing is in 
the path of such probability streams, it may have a very 
large probability, even if the energy barrier eq associated 
with this packing is low. 

We note that in our intuitive picture, probability flux 
is akin to water flow in a rugged mountainous landscape. 
This analogy, however, may be oversimplified due to the 
highly multidimensional character of the PEL for partic- 
ulate systems. 



B. Displacement fluctuations 

1. Displacement and dynamical matrices 

After a digression on possible non-local features of the 
PEL that can control the large fluctuations in the prob- 
abilities of MS packings, we now return to our analysis 
of local features of the PEL near a given MS packing. 
In this section we focus on measuring the shape and size 
of the local region visited by the system when the MS 
packing is thermally excited. 

In general, the local shape of the PEL basin can be 
studied using both static and dynamic techniques. The 
static method relies on an analysis of the eigenvalue spec- 
trum of the dynamical (Hessian) matrix|3|and0] The dy- 
namic technique that is applied in this section is based 
on the evaluation of the 2A'^-dimensional matrix of dis- 
placements away from the reference MS packing while 
the system fluctuates after input of thermal energy. The 
displacement matrix is deflned as 

AajT? = {{ria - ^ao) {rjf} - r-jpo)), (12) 

where riao are the coordinates of the reference MS pack- 
ing, and Via are the current coordinates. Again, the in- 
dices i and j refer to the particles, and a, (3 = x,y rep- 
resent the Cartesian coordinates. The average in Eq. IT?! 
is taken both over times t < t less than the minimal 
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FIG. 14: Comparison of the eigenvalues of the displacement 
matrix di/T for several temperatures T (curves) and those of 
the inverse dynamical matrix l/irii (open squares) for (a) a 
system that has been compressed by 1% from the jamming 
threshold and (b) a MS packing. Note that when plotting the 
eigenvalues of the dynamical and displacement matrices, we 
include only the dN — d nontrivial eigenvalues and order them 
from smallest to largest. 

breaking time and over at least 20 realizations at each T, 
weighted by the corresponding breaking time r. 

In harmonic systems, the eigenvalues of the displace- 
ment matrix di are trivially related to those of the dy- 
namical matrix 

= — . (13) 

Near an MS packing, however, our system of disks is 
strongly anharmonic, due to the one-sided character of 
the spring potential at the jamming threshold. Thus, 
relation 1131 is not guaranteed to hold even in the low- 
temperature limit. 

We find that for compressed systems with finite parti- 
cle overlaps there is a characteristic temperature Tc be- 
low which relation 1131 is satisfied. For example, a sys- 
tem compressed by 1 % above the jamming threshold be- 
haves harmonically for T < Tc ~ 10~^, as illustrated in 
Fig. EI a). However, the temperature Tc tends to zero as 
the amount of overlap decreases; therefore relation [T51 is 
not valid for MS packings even at infinitesimal tempera- 
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FIG. 15: (a) The minimum dmin, (b) median dmcd, and (c) 
maximum dmax eigenvalues of the displacement matrix plot- 
ted versus the minimal energy-barrier height eo at temper- 
ature T/eo = 0.2 for A*' = 10 and all 105 packings in our 
sample. The eigenvalues of the displacement matrix scale as 
di ~ Eg' , and the power-law exponents 7i are given in Fig. Ilfcil 
For dmax, a few packings in the subset {Pfc}''^'^'^ (solid circles) 
deviate from the main trend. 



tures. The quantities on the left and right hand sides of 
Eg. 1131 can differ by more than an order of magnitude for 
MS packings as shown in Fig. I14r &). 

In what foUows, we focus on the displacement matrix 
1121 (rather than the dynamical matrix) because it more 
reliably measures the local shape of potential energy 
landscape near MS packings. Moreover, at finite tem- 
peratures the matrix Diajp yields information about the 
entire region of configuration space with energy V < eQ 
near a given MS packing. Note that since we are consid- 
ering thermally fluctuating systems in which all particles 
move, the distinction between rattler and non-rattler par- 
ticles becomes less important and we will not distinguish 
between these two types of particles in the following dis- 
cussion. 
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FIG. 16: The power-law exponents 7i that determine how the 
eigenvalues of the displacement matrix scale with the minimal 
energy-barrier height eo. The index i runs from the smallest 
eigenvalue of the displacement matrix to the largest. For har- 
monic systems, all 7i = 1. 



2. Eigenvalues of displacement matrix 

The results displayed in Fig. E| (b) show that the 
eigenvalue spectrum of the displacement matrix for MS 
packings can be highly nonuniform. In this example, the 
ratio of the largest to smallest displacement eigenvalues 
rfmax/rfmin IS approximately 10"^ at low temperatures. It 
is thus interesting to examine how this ratio, and more 
generally the shape of the basin near each MS packing, 
varies from one packing to another. 

The shape of the local basins was explored by evaluat- 
ing the displacement matrix^jfor different MS packings 
at a fixed value of T/eo, so that for each packing we 
study comparable relative displacements away from po- 
tential energy minimum. We find that the eigenvalues of 
the displacement matrix exhibit power-law scaling with 
the energy-barrier height. 



■=0 ' 



(14) 



where the scaling exponent 7i depends on the position i 
in the set of d{N — 1) eigenvalues ordered by their mag- 
nitude. 

The power-law behavior^] is illustrated in Fig. 1151 for 
the minimum, median, and maximum eigenvalues of the 
displacement matrix for N = 10 at fixed T/eo = 0.2. We 
also performed measurements at T/eo = 0.3 and obtained 
similar results. The power-laws extend over at least 10 
orders of magnitude in eo , although there is some scatter 
in the data for the largest eigenvalue dmax- c'max for 
several of the packings from the most frequent subset 
{-Pfe}«=°-^ deviate fr om the main trend. 

The scaling exponents are given in Fig. 1151 and range 
from approximately 0.6 for the largest eigenvalue to near 
1.0 for eigenvalues below the median. The fact that the 
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FIG. 17: Ratio dmax/rfmin of the maximum and minimum 
eigenvalues of displacement matrix versus the minimal energy- 
barrier height eo at temperature T/eo = 0.2 for A'^ = 10 and 
all 105 packings in our sample. The long-dashed line has slope 
—0.5. The most frequent packings in the subset {Vk}''^'^'^ are 
indicated by filled circles while the rest are indicated by open 
circles. 



largest several scaling exponents differ significantly from 
unity again indicates that MS packings are anharmonic, 
since in harmonic systems ji — 1 (assuming that the 
eigenvalues of the dynamical matrix are independent of 

To illustrate the non-uniformity of the displacement- 
matrix spectrum, in Fig. 1171 we plot the ratio rfmax/c'min 
versus eq at fixed T/eq = 0.2. We find that the ratio 
roughly obeys the power-law 



^max -0.5 



(15) 



consistent with the results shown in Figs. El (e^) and 
(c) and El Fig. El emphasizes that MS packings can 
possess extremely nonuniform displacement fluctuations 
with draax/ drain as large as 10^ for packings with small 
energy barriers. As in Fig. El (c), there are several fre- 
quent MS packings with eo < 10^^ that have consider- 
ably larger ratios than this trend, but the vast majority 
of points obey (fT5|) . 

Our data shows that MS packings possess highly 
nonuniform displacement eigenvalue spectra (nearly all 
have (imax/<^min > 100) and the non-uniformity sub- 
stantially increases with decreasing eo- This suggests 
that the basins associated with MS packings are highly 
anisotropic in configuration space, and this gives rise to 
displacement fluctuations that are much larger in one or 
a few directions than others. Thus, one simple picture 
is that MS packings break only along particular direc- 
tions in configuration space and that the breaking direc- 
tion will be correlated with the direction in which the 
displacement fluctuations are the largest. We will inves- 
tigate this intuitive picture in future studies. 
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FIG. 18: Definition of contact angles. 6jk is the angle be- 
tween the lines that connect the central particle i to a pair 
of adjacent, contacting neighbors j and k. There are four 
contact angles associated with particle i. 

The results presented in Figs. I15H17I allow us to deter- 
mine if there is a direct link between the volume Qk and 
the packing probabilities fk, where D,k is the volume of 
configuration space near a given MS packing that con- 
tains points with potential energy y < eo . An estimate 
of this volume can be obtained by assuming that at tem- 
perature T ~ eo the system explores a large portion of 
Accordingly, we find 



d{N-l) 
i=l 



(16) 



where r] — ^ J2iLi 7i- The last expression in was 
obtained using the relation between di and eo in (|14|l . 

Thus, if we assume that the packing probabilities are 
controlled by the volume Q^, we obtain 



fk 



(17) 



where 77 « 7 for the 10-particle system considered here. 
However, our results in Fig.ll2lshow that the dependence 
of the packing probabilities on ep is much weaker than 
that predicted by 1)17(1 . The exponent for the 10-particle 
data is approximately 2, not 7. 

Thus, our results indicate that the packing probabili- 
ties are determined by a much lower dimensional quan- 
tity than the volume flk ■ Our results are consistent with 
a scenario in which the packing frequencies are deter- 
mined by only the largest several displacement eigenval- 
ues. It is likely that the packing probabilities are corre- 
lated with the largest displacement eigenvalues because 
they may correspond to the (initial) escape directions 
from the basin of a MS packing. 



V. STRUCTURAL PROPERTIES OF MS 
PACKINGS 

We also examined structural properties of MS packings 
in an attempt to identify features that control MS pack- 
ing probabilities. Specifically, we measured the proba- 
bility distributions for the angles between lines connect- 
ing centers of particles in contact (cf., the definition in 
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FIG. 19: (a) Maximum contact angle ^max (in degrees) for 
each MS packing versus the frequency with which it oc- 
curs for a 10-particle system using the MD energy minimiza- 
tion method, (b) The distribution of the largest contact angle 
^max (in degrees) on each particle for the 100 most frequent 
(thin lines) and 100 most infrequent (thick lines) for the same 
system in (a). The dashed (solid) lines include all (larger half) 
^max in each configuration. 

Fig. I18|) . We focused on the distribution of maximal an- 
gles between particle contacts because large angles (close 
to 180°) correspond to unsupported (and thus unstable) 
chains of nearly collinear particles. 

As defined in Fig. 1181 9jk is the angle between the lines 
connecting the central particle i to pairs of adjacent, con- 
tacting neighbors j and k. A given particle will possess 
ric contact angles, where Uc is the number of contacts for 
that particle. We calculated the maximum contact an- 
gle ^Jnax fo'" each of the N' non-rattler particles and the 
maximum angle ^max for each configuration. 

In Fig.^l(a), we present a scatter plot of 6'niax f^or each 
MS packing in a system with = 10 particles versus the 
frequency fk with which the packing occurs for the MD 
packing-generation protocol. This plot shows that there 
is no strong correlation between B^ax and fk- In fact, 
some highly probable MS packings possess large contact 
angles near 180°. 

Unsupported nearly linear chains of particles imply a 
low energy barrier. Since the barrier heights are corre- 
lated with the probability fk for the MD protocol, our 
results may thus suggest an important role of near con- 
tacts between particles. Such contacts do not directly 



support the chain but they may prevent transitions to 
other MS packings when the system is fluctuating [47l| . 
We will return to this problem in future investigations. 

We also measured the distribution of the largest con- 
tact angles on each particle ^^ax- Fig. ^1 (b), we 
compare the distributions of ^^ax for the 100 most fre- 
quent (thin lines) and 100 most infrequent (thick lines) 
MS packings generated for a 10-particle system using the 
MD protocol. The figure shows two sets of curves. The 
probability distribution represented by the solid lines in- 
cludes all ^Jj^ax ill the system; the dashed hues represent 
the distribution of only the largest half of the angles for 
each MS packing. 

In contrast to the results represented in Fig. ^] (a), 
the distributions depicted in Fig. ^| (b) do show a clear 
correlation between large angles ^^ax and the frequencies 
fk of MS packings: The infrequent ones have an excess 
of large angles near 160° compared to the frequent pack- 
ings. Thus infrequent MS packings tend to have multiple 
particles with large contact angles. 

The structural differences between the frequent and 
infrequent states that we see in Fig. E| are quite sub- 
tle. We do not have any significant correlation with the 
single angle ^max in the packing. The differences show 
up, however, collectively in the distribution of the con- 
tact angles for each particle in the packing. One of our 
long-term goals is to connect important geometrical fea- 
tures of configuration space to structural properties of 
MS packings that can be measured in experiments. 



VI. IS THERE A HIDDEN RANDOM 
VARIABLE? 

In previous sections we presented a detailed study of 
the probability distribution of MS packings in small 2d 
systems of bidisperse frictionless disks. We showed that 
the probabilities of individual packings may differ by or- 
ders of magnitudes, not only as function of packing frac- 
tion, but also in individual narrow packing-fraction inter- 
vals. Moreover, we identified important features of the 
packing probabilities fk that are only weakly affected by 
the details of the packing-generation protocol. 

Since our packing-generation algorithms do not target 
any specific packings, there must exist important prop- 
erties of the multidimensional configuration space that 
give rise to such widely varying packing probabilities. We 
have examined several local properties of the PEL near 
the MS packing configurations and have found gross cor- 
relations between these properties and the packing prob- 
abilities. However, there is large statistical scatter in the 
data, and thus we conclude that none of the local fea- 
tures of PEL that wc have examined can fully explain 
the packing probability distribution. 

Yet, in spite of the complexity of the problem we have 
found several interesting regularities. In particular we 
have determined that the probabilities fk in a narrow 
packing-fraction interval, evaluated for different values 
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of N, can be rescaled onto a single master curve with the 
characteristic shape shown in Fig. |Sl(b). In this section 
we further explore this striking similarity of the sorted 
probability distributions. We base our analysis on a sim- 
ple phenomenological model that has been inspired, in 
part, by the correlation between the packing probabili- 
ties and the minimal energy barriers [cf. Fig. 1121 (h)]. 

In this model we characterize each MS packing by a 
set of M independent continuous random variables 



Xi,...,XM > 0, 



(18) 



all with the same probability distribution tti (x) , which is 
approximately uniform for x < Xmax and quickly decays 
beyond Xmax- The number of random variables Xk is 
comparable to the number of degrees of freedom in the 
system M = 0{Nd). Our central assumption is that 
the probability fk of a MS packing within the packing- 
fraction interval of interest is controlled by the smallest 
of the random variables 1181 



fk min(a;i, . . . ,xm)- 

Using this assumption, the probability density tt 
the random variable 



y = min(xi, . . . ,xai) 



(19) 

mm for 

(20) 



can be written as 



where 



^min(2/) = Af7ri(y)nf-i(2/) 

- nf (,), 

ay 



Ili{y) = / TTi{x)dx 

J y 



(21) 



(22) 



is the probability that x > y. 

From this result we can evaluate the expected value of 
the number of MS packings k{y) that have < y < y, 



TTmUy)dy = i-^f{y). 



(23) 



where fcmax is a total number of MS packings. 

Now we apply the above results to the packing fraction 
interval in which we have fcmax = f^s(A(/)) states. After 
inserting assumption ll9l into l23l wc obtain the expression 



i-nf(a-iA), 



(24) 



where a is the proportionality constant in Eq. 1191 Ac- 
cordingly, for a given probability distribution Hi, our 
analysis yields the relation between the sorted probabil- 
ities /fc and the index k in the sorted sequence of states. 
This relation corresponds to the plot shown in Fig. El 
(strictly speaking, to its inverse). 
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FIG. 20: The right-hand side of Eq. (1 - fc/fcmax)^''", 
with M — dN — d, is plotted versus the normalized discrete 
packing frequencies fk / /™'"^ over the narrow packing fraction 
intervals in Fig. |5] using MS packings from the CG method 
in A'' = 10 (dotted line), 12 (long-dashed line), and 14 (solid 
line) particle systems. 



Our theory cannot be directly verified without specify- 
ing the distribution of the random variables Xk, and we 
do not have any a priori information regarding this dis- 
tribution. Our goal here, therefore, is more limited. We 
simply want to determine whether or not the numerical 
results shown in Fig. |^ are consistent with our assump- 
tion that 7ri(a;) is a relatively uniform function in some 
range of x, outside of which it quickly decays. 

In making this assumption, we have in mind features of 
the PEL such as the distance from a given MS packing to 
the passes in the rim of the local potential-energy basin. 
In each direction, the distance to the rim is smaller than 
the particle diameter; it is also conceivable that the dis- 
tance to the closest pass determines the probability fk- 
However, the specific meaning of the hypothetical ran- 
dom variables ^1 in our phenomenological theory at this 
point has not been fleshed out. 

To determine the approximate form of the probability 
distribution Hi from our numerical results, relation 1241 is 
inverted. 



ni(a-Vfc) = (1-fc/fcn 



sl/M 



(25) 



and the data sets represented in Fig.|H|are replotted in a 
form that emphasizes the structure of our model. Specifi- 
cally, in Fig. 1201 the quantity on the right-hand side of Eq. 
I2S1 is shown versus the normalized probability fk / fk^'^^ 
48]. The results in Fig. [201 indicate that the transformed 
quantity 1211 is a nearly linear function of fk- By Eos. 1221 
and 1251 this linear behavior is consistent with 



7ri(x) ^ <d{x), 



(26) 



where 8(a;) is the Heaviside step function. The form 
of the probability distribution tti determined from our 
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numerical results is thus compatible with assumptions of 
our model. 

The results in Fig. [501 cannot be treated as a direct 
verification of our model, and the model itself is rather 
ad hoc. Nevertheless, the simplicity of the result 1211 is 
notable. We thus believe that our model captures at 
least some essential aspects of the problem. However, 
the exact source of the very broad distribution of the 
probabilities fk and the nature of the self-similarity of 
this distribution for different values of N (as revealed by 
the results shown in Fig. ^ still remains an important 
open problem. 



VII. CONCLUSIONS AND FUTURE 
DIRECTIONS 

We have performed extensive numerical simulations 
with the aim of generating mechanically stable (MS) 
packings of frictionless disks in small bidisperse systems. 
The MS packings are created using a protocol in which we 
successively grow and shrink soft, purely repulsive parti- 
cles. Each compression or decompression step is followed 
by potential energy minimization, until all particle over- 
laps are infinitesimal. We focus on small systems with 
at most 14 particles in 2D because in these systems we 
are able to find nearly all distinct MS packings and can 
therefore accurately measure the frequency with which 
each packing occurs. 

One of the principal results in this work is that MS 
packing frequencies differ by many orders of magnitude 
both as a function of packing fraction and within narrow 
packing-fraction intervals. We have implemented here a 
fairly generic algorithm for generating MS packings, and 
this algorithm does not specifically target any of them. 
Yet we find that packing frequencies are extremely varied; 
moreover, the frequency variation increases with system 
size. We also find that the probability distribution can 
be scaled onto a single master curve. 

We argue that these striking results are important in 
a broader context of theories of jammed granular me- 
dia and glassy materials. In thermodynamic theories 
for dense granular media [s^ it is usually assumed that 
MS packings within a small packing-fraction interval are 
equally probable. For our packing-generation protocol, 
this assumption is certainly not valid. Although we do 
not yet have sufficient data, we believe that MS packings 
will also not be equally probable for other commonly used 
experimental protocols, e.g., slow shear or vibration 
under gravity ^5Qj. However, thermodynamic theories 
based on the equal-probability assumption are often ap- 
plied to understand the properties of sheared or vibrated 
granular materials 34,51]. 

Our findings show that the often-used equal- 
probability assumption for stable grain packings in gran- 
ular matter and inherent structures in glass-forming liq- 
uids should be re-examined. If it turns out that the as- 
sumption is generally violated, except for some unphys- 



ical, highly specialized algorithms (and we expect that 
this is the case), thermodynamic theories of disordered 
granular packings will need to be significantly reformu- 
lated. 

In this work we focused entirely on a system of fric- 
tionless particles. However, we would like to point out 
that for frictional particles there is another conceptual 
difficulty with the assumption of equal-probability pack- 
ings. Since static friction can arrest particle motion at 
different contact angles (analogous to a block that can 
stop at any position on a wedge) MS packings of fric- 
tional particles do not form points in configuration space, 
but rather continuous hyper-surfaces. Since packings of 
frictional particles are often hyperstatic |53 |. the dimen- 
sionality of these hyper-surfaces changes from packing to 
packing, and it is thus difficult to introduce an appropri- 
ate probability measure. 

Returning to the summary of the key results of our 
study, we note that important features of the MS packing 
probabilities do not change when we alter the packing- 
generation protocol. Thus, we argue that protocol- 
independent properties of configuration space must play 
an important role in determining the MS packing fre- 
quencies. To investigate the connection between geo- 
metric properties of configuration space and MS packing 
probabilities, we added thermal energy to a set of MS 
packings and then measured several quantities as the sys- 
tem fluctuated. We monitored the time that elapsed be- 
fore a MS packing broke to a distinct one, which allowed 
us to determine the heights of energy barriers that sepa- 
rate one MS packing from another. We also studied the 
displacement fluctuations in all possible directions away 
from the original MS packing to infer crucial features of 
the shape of the basin near each packing. 

We found a gross correlation between the frequencies 
fk of MS packings and the height of the lowest energy 
barrier eq separating this packing from other ones. The 
MD packing frequencies roughly scale as fk ~ £q with 
A ~ 2, but there is significant scatter in the data. In 
addition, we found that the eigenvalues of the displace- 
ment matrix scale as di ^ eJJ' with 0.6 ^ 7i ^ 1- These 
results suggest that the MS packing frequencies are de- 
termined by one or a few degrees of freedom, not by the 
local volume of configuration space near each basin (in 
which case fk would scale much more strongly with eq). 
However, the scatter in our data implies that there are 
important unknown variables that are linked to the MS 
packing probabilities. 

Our results clearly indicate that this is complex prob- 
lem and much more work needs to be done to fully under- 
stand what determines the MS packing frequencies. Here 
we mention briefly some directions that we are actively 
pursuing to address this question. 1) We are measuring 
the hypervolumes of the regions in configuration space 
whose vertexes are the nearby low-order saddle points of 
each MS packing. We will investigate the relation be- 
tween these hypervolumes and packing probabilities. 2) 
We are studying the {dN — c?)-dimensional breaking vec- 
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tor that connects the initial MS packing to the MS pack- 
ing to which it breaks. We want to determine whether 
or not the breaking vector is correlated with the direc- 
tions of large displacements when the system is thermally 
fluctuating. 
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